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| A quantum master equation (QME) is derived for the many-body density matrix of an open 

. current-carrying system weakly coupled to two metal leads. The dynamics and the steady-state 

X^/y ' properties of the system for arbitrary bias are studied using projection operator techniques, which 

keep track of number of electrons in the system. We show that coherences between system states 
with different number of electrons, n, (Fock space coherences) do not contribute to the transport 
to second order in system-lead coupling. However, coherences between states with the same n 
may effect transport properties when the damping rate is of the order or faster then the system 
Bohr frequencies. For large bias, when all the system many-body states lie between the chemical 

■ potentials of the two leads, we recover previous results. In the rotating wave approximation (when 
the damping is slow compared to the Bohr frequencies of the system) , the dynamics of populations 
and the coherences in the system eigenbasis are decoupled. The QME then reduces to a birth and 
death master equation for populations. 
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I . I. INTRODUCTION 

Electron transport through a quantum dot (QD) or a single molecule has recently received considerable experimental 
0- 13- E3- Q- 0- E3] and theoretical 0, 0, fToL fllL \l$L Il3j attention. The progress in the fabrication of devices such as 
quantum dots (whose size and geometry can be controlled with high precision |14|) or single molecule junctions, makes 
it possible to investigate quantum effects on transport and provides a test for methods of nonequilibrium statistical 
, mechanics. In analogy with laser optical spectroscopy |15t Il6| - electron transport through a quantum system (QD or 
single molecule) can be used to p robe its nonequilibrium properties through the current-vo ltage (I IV) characteristics. 
The scattering matrix (SM) |l7llTs| and the non-equilibrium Greens function (NEGF) 0, 12fll21| methods have been 
used for predicting I/V characteristics of quantum systems connected to two metal leads. Both theories are exact 
in their respective domains: SM is limited to elastic processes while the NEGF can treat both elastic and inelastic 
processes. 

The quantum master equation (QME) approach is an alternative tool for studying the irreversible dynamics of 
quantum systems coupled to a macroscopic environment [U |2jJ . Owing to its simple structure, it provides an 
intuitive understanding of the system dynamics and has been used in various fields such as quantum optics |25| , 
solid state physics |2q , and chemical dynami cs 1 1 5lL Recently, it has been applied to study electron tunneling through 
molecules or coupled quantum dots j2j, IH, I29L l3Ct l3ll l32| . Fransson and Rasander [33J have recently used a QME 
approach to study the rectification properties of a system of coupled QDs by analyzing the occupation of two-electron 
triplet states as a function of the ratio of the interdot coupling and the energy splitting between the two QDs. 

Gurvitz and Prager were the first to derive a hierarchy of QMEs which keeps track of the number of electrons 
transferred from the source-lead to the collector-lead. Using this hierarchy they studied the effects of quantum 
coherence and Coulomb blockade on steady state electron transport in the high bias limit. In this limit all energy 
levels of the system are below the chemical potential of the left lead (source) and above the right lead (collector) (Fig. 
I ). The relevant Fermi functions for the left and the right leads are then unity and zero, respectively. Rammer et al 
29] have used a QME to describe the direct tunneling (where the system never gets charged) in quantum junctions. 
Recently Pedersen and Wacker |3(j have generalized the standard rate equation and included approximately two- 
electron transfer processes by going beyond the second order perturbation in system-lead coupling. 
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In this paper, we use projection operator techniques to derive a new hierarchy of QME for the many-body density 
matrices p n representing the system with n electrons. Electron transport through a quantum system is expressed 
in terms of the four processes describing the charging (a + and b + ) or discharging (a~ and b~) of the system at the 
left and the right leads (Fig. 1). Yan et at ,3l| have used the same model to compute the steady-state current in 
the system by keeping track of the number of electrons at the collector-lead. In the limit of large bias, when the 
backward transport (corresponding to electron moving in the direction unfavored by the potential difference between 
the leads) can be neglected, we recover their results. Otherwise, it is necessary to identify n as the number of electrons 
in the molecule, as done here. We solve our equations for a model system of two coupled quantum dots and study 
the effect of quantum coherences on electron transport. Coherence effects in quantum junctions have been studied 
in the past. Using the scattering matrix approach, Sautet et al |34j have found interference effects on the scanning 
tunneling microscopy images of molecules adsorbed on a metal surface. These effects arise from coherences between 
different open channels for the tunneling electrons. 

The paper is organized as follows: In section^] we present the Hamiltonian and define a projection operator which 
keeps track of the system's charge. We derive the QME and discuss its connection with earlier works. In section ITTll 
we show that under different approximations our QME recovers previous results and assumes a Lindblad form. Under 
the rotating wave approximation (in the system eigenbasis) the QME provides a very simple single particle picture 
of the dynamics of populations and coherences. In section llVl we study the dynamics, current and the charge of two 
coupled quantum dots (QD). In section we present numerical results and discuss the effect of quantum coherences 
on the current. Conclusions are drawn is section IVTl 

II. THE QUANTUM MASTER EQUATION 

The Hamiltonian of a quantum junction is given by the sum of the Hamiltonians for the isolated system, H s , the 
left, Hl, and the right leads, Hr, and the lead-system coupling (Ht)- 
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= H S + H L + H R + H T 
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Ht 


= ^ [ T s" c i c " + h.c] 

sv 


(5) 



where s, I and r represent system, left, and right lead orbitals, respectively, and v = I, r. T s i and T sr are the transfer 
coupling elements between the leads and the system. Direct coupling (tunneling) between the leads is neglected 29] . 
c\(c s ), cj(ci) and <4(c r ) are electron creation (annihilation) operators which satisfy the Fermi commutation relations 

{c k ,cl,} = 6 k k>, {4>4'} = i c k,Ck'} = 0, k,k' = s,l,r (6) 

where {A, B} = AB + BA. 

The many-electron eigenstates of the system form a ladder of manifolds: the n'th manifold \np) contains the states 
(p) with n-electrons. Each interaction with the leads, Eq. ©, can change n to n± 1. The total many-electron density 
matrix in Fock space can be expanded as, 

Pt= £ (Cn\np){mq\. (7) 

n,m,p,q 

The diagonal (n = m) blocks represent Fock space populations (FSP) of system states with n electrons whereas the 
n m blocks are Fock space coherences (FSC). When the system is brought into contact with the leads, it is initially 
in the n'th FSP block, p^ n . As time evolves, FSC are developed, inducing transitions to other FSP blocks nil, 
n ± 2, etc. At steady state, these blocks reach a stationary distribution and the current through the system can be 
calculated using time derivatives of the FSP. 

Our first step is to derive a quantum master equation in Fock space which keeps track of the FSP and eliminates all 
FSC by incorporating them through relaxation kernels 2m- The Markovian master equation holds when the dephasing 
rates of FSP are large. In that case the steady state coherences are small, and progressively decrease for higher order 
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coherences, i.e. as \n — m| in Eq. Q is increased. The dominant terms are m = n ± 1 and the master equation rates 
can be calculated to second order in the coupling (T) with the leads. As the FSC dephasing rates decrease, one should 
calculate the rates to higher order in T. This problem is formally equivalent to the multiphoton excitation of molecules 
or atoms; the molecular states are divided into n-photon manifolds, and n-quantum coherences are eliminated to derive 
a Pauli master equation for the populations. Coupling with the radiation field in the rotating wave approximation 
plays the role of the coupling with the leads. The time-convolutionless projection operator techniques developed for 
multiphoton processes |35j| can be applied towards the calculation of molecular currents. 

To derive a reduced description in the system space, we define the projection operators V n which act on the 
many-body wave function "J 29] 

V n ^{r 1 ,---r N )=e n (r 1 ,---r N )^(r 1 ,---r N ) , (8) 

where O n (ri, ■ ■ - rjv) = 1 if precisely n space-points belong to the system subspace and it vanishes otherwise. V n is 
thus a Fock space projection operator onto states with n electrons in the system. The projected many-body density 
matrix (p n ) onto the system subspace with n electrons is defined as, 

p n = Tv lead {V nPT V n }. (9) 

Note that by defining the projection operator for a fixed number of electrons in the system, we ignore the coherences 
between the leads and the system. The projection of the many-body density matrix can also be formulated in Liouville 
space using the projection superoperator (C„) as C n px — p n PB, where p B is the density matrix of the leads (bath). 

The QMEs for the projected many-body density matrices of the system is derived in Appendix A using second 
order perturbation theory in the system-lead coupling and by treating the two leads as infinite electron reservoirs: 

= -i[H s ,p n (t)} + ]T [a ss ,c sl p n+1 (t)4-p ss ,p n (t)c sl ci 

ss' 

- a ss ,clc s ,p n (t)+(3 ss ,cip n ~ 1 (t)c s , +h.c] , (10) 

with 

T SV T*, V {1 - f v ) 



OO 



dre^' T a ss , (r) = lim V sv 8 ' vK Jv> (11) 



/ dre^'^(r) = lim V — 



L, I 'It, - ' r .l-iD : lim > ] T^hk , (12) 

' - . + IT] 



where // (f r ) is the Fermi distribution of the left (right) lead with chemical potential pl = Po + (/J-r = Mo) and 
eV is the applied bias. a SS '(r)'s and /3 ss '(r)'s are the lead correlation functions, Eq. I|A16I) . and have the symmetry 

a* ss ,(r) = a s , s (-r) , &V(r) = M~r) ■ (13) 
Taking into account that the leads are macroscopic and have a continuous density of states, Eq. (|A16|I gives 

^2e-^T su T;, v f u - Wde n x {e)e— (e)T^ X > (e)f x (e) , (14) 

v X J 

where X = L,R and nx is the density of states of lead X. Substituting the relation 

dre lUT =tv8(u)+iV- , (15) 



(I 



in Eqs. (|llfl we obtain 

a ss ' = ; fa=Y,fi? ( 16 ) 

x x 



<*& = ™,( es ,)rf)(e,)Tf'*( £ ,)(l-/x( £ .)) +I / ^ n ^ X) ^^-f^) (i 7 ) 

(18) 



a {X) , srr(x) f \rr(X)*, s, f s. f , nx{e)T ( s X \e)TlP*{e)f x {e) 
Pi/ = 7rn x (e s /)Ti >(es')T^ {e s ,)f x (e s >) + i / deV 
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The real parts of a sa > and f3 ss i define the system to leads and leads to system electron transfer rates, respectively. The 

i x\ 

imaginary parts represent level shifts. ^R.a ss is the rate with which electrons are transferred from the s'th system 

( x ) 

orbital to lead X while $t/3ss is the transfer rate of electrons from lead X to the system orbital. Thus $la ss > and 
3?/3 ss ' are associated with the processes where the system undergoes transition between many-body states which differ 
by a single electron. When the external bias is large enough so that /i(e s ) = 1 and fn(e s ) = 0, Sfta^, = $lf3 R , = 0. 
This means that the backward flow of electrons (electrons moving against the applied bias) from the right to the 
left lead vanishes and that each time the number of electrons in the system decreases it corresponds to an electron 
tunneling to the right lead. Keeping track of the electrons in the system is therefore directly related to counting of the 
electrons collected in the right lead. In such case, we recover the result of Yan jsj. However, in general it is essential 
to recognize that p n is the density matrix of the system with n electrons residing in the system. When n decreases 
by one, the electron will, with a higher probability, be collected in the right lead, but could also be collected in the 
left lead. The effect of this last process on the dynamics is not captured in the other QME [3lJ but is made clear in 
our QME by the use of projection operators. 

To appreciate the reduction involved in the QME, let us consider a system with n electrons and M(n < M) orbitals. 
The number of n-electron many-body states is then C^f = ntjj^rcj and the total number of many body states is 

JVtot(M) = J2n=0 C n = 2 M Thc fml many-body density matrix is N to t(M) x N to t(M). Because the FSC between 
many-body states with different n are eliminated, the size of the reduced many-body density matrix is Xl^Lo^n^) 2 - 
In the full Liouville space of the system, the many-body density matrix is an 7V t 2 ot (M) vector. However, the projected 

many-body density matrix p s = J2 n P n m this space contains (J2 n =o ^n) 2 ~ J2 n =o(^n) 2 elements which are zero. 
Our QME is therefore defined in a reduced many-body Liouville space of the system where the FSC have been 
eliminated. 



III. LIMITING CASES AND THE LINDBLAD FORM 



A. High bias limit 



When all many-body states of the system lie within the chemical potentials of two leads, the reverse electron 
tunneling can be ignored since the Fermi functions for the left and right leads are /z,(e s ) = 1 and fnX^s) = |28l l36| . 
If we neglect the principal parts in Eqs. (|17|l and (|18|) . the matrices a and (3 become hermitian. In this case, since 



<> = Pfs> 

Ojs' and f3 ss . 
[37| . we get 



we have a s 
Defining B s = 



= n R T R T R and /3 SS , = n L T^. We have further ignored the energy dependence of 
^T R c s: D s = v / n r T/'*c| and summing the QME (JTUJ over n so that p(t) = £ n p n (t) 



= -i[H s , p]+J2 [ 2B s'PB\ - B\B s ,p - pB\,B s + 2D s ,pD\ - D\D s ,p- p{t)D\,D s 



(19) 



Keeping in mind that FSC are zero so that terms, such as B s >pD\, B\,pB s , etc. which lead to FSC, vanish, Eq. ijT§|) 
assumes a Lindblad form, 



p = -i 



i[H s , P } + J2 [zAs>pAi - A\A s ,p - P A\,A S 



(20) 



where A s — B s + D s = V n R T^c s + vn^T^ c\. Gurvitz [2j| has studied the effect of coherences in a QD system 
connected in series in the high bias limit. In Appendix B we show that Our QME reduces to Gurvitz equations in 
this limit. 



B. The Rotating Wave Approximation 

When the interaction between the system and the leads is weak enough for the damping effects to be slow compared 
to the Bohr frequencies of the system, we can simplify Eq. (|10Jl by using the rotating wave approximation 0, |22, H3| . 
This approximation is often performed on the Markovian form of the Redfield equation in order to prevent the 
possible breakdown of positivity [2^, 123. 138L |3!| due to the nonMarkovian effects of thc initial dynamics |4Ct l4ll |42| . 
Transforming the master equation to the interaction picture, we get Eq. <|A15ll with the Markovian approximation 
J * dr — > J °° dr and with the Born approximation p n (t — r) — > p n (t) . Since the damping rate of the density matrix 
in the interaction representation is small compared to the Bohr frequencies of the system, we can time average 
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(limT->oo j! Jq dt) the fast oscillations due to the terms e ie "'* in Eq. (|A15|) . This allows to eliminate the nondiagonal 
elements of the correlation function matrices [a ss > — a ss 5 ss i and ss i — [3 SS 5 SS ']. Going back to the Schrodinger 
picture, our equation reads 



P 



■i[H s , p n ] + [a ss c sP n+1 ct - & ss p n c s c\ - a ss c) s c s p n + p ss c\p n - l c s + h.c] . (21) 



By projecting this equation into the reduced Fock basis, we find that the coherences are decoupled from the popula- 
tions. 

When the level shifts in Eqs. (|17|l and IjlSjl are ignored, the matrices a and (3 in Eq. Ij21(l become Hcrmitian. Using 
a s = y/a S sC Sl b s — \J (3sscl and following the same steps as in case of high bias limit, we find that Eq. (|21|) when 
summed over n is of a Lindblad form similar to Eq. I|2UI) : 

P = -i{H s ,p}+J2[2A s pAl-AlA s p-pAlA s ] . (22) 

s 

Note that in the RWA, since the matrices a and f3 are diagonal, the sum in Eq. Il22[l only runs over one index s and, 

L(R) 

unlike the high bias case, the coupling coefficients T ss y can be energy dependent. 

Thus, we conclude that both in the high bias limit and the RWA form our QME are of Lindblad form which guarantees 
to preserve the positivity and the hermiticity of the density matrix. 

The dynamics of the reduced many-body density matrix [Eq. (|21|l ] can be expressed in terms of the time evolution 
of the single-orbital density matrix p s corresponding to the sth orbital. 

Ps = -ie ss [clc s ,p s ] + 2$ta ss c sPs cl + 2$tf3 ss clp s c s 

-a ss c\c s p s - a* ss p s c\c s - f3* s c s c\p s - (3 ss p s c s c\ , 

This means that if the initial many-body density matrix in Fock space is a direct product of single-orbital density 
matrices p = pi ® p2 <8> • • • <S> Pm it will remain so at all times. However, even if the initial many-body density matrix 
is not a tensor product, since in the RWA, the r.h.s. of Eq. I|21|) is a sum of contributions from various orbitals, we 
can still factorize the many-body populations as products of single-orbital populations 

AI 

(»!•■• n M \p\ni ■ ■ ■ n M ) = J[ p# (23) 

s=l 

where p^f 1 — (n s \p s \n s ) and n s is the occupation of sth orbital. The dynamics of the single-orbital occupation is given 
by 



a,- ; 2 [ -^wLH ' (24) 



We can readily find the steady state distribution 



W = MPss (s) = «Q M , , 



The many-body steady state distribution can be directly obtained using (PI) and {25J. 



IV. MODEL CALCULATIONS 



We consider a system of two coupled quantum dots (QD) each having a single orbital in the energy range between 
the chemical potentials of the left (pl) and the right (pr) leads. Depending on the interdot coupling, the system 
orbitals may either be localized (weak coupling) or delocalized (strong coupling). The system Hamiltonian is 

H s = eicjci + e 2 4c 2 , (26) 

where t\ and e 2 are system orbital energies and we have ignored the charging effects due to electron-electron inter- 
actions (Coulomb-blockade) [28L l43l| in the system. We denote the many-body states {|ni,ri2)}, where n% and n 2 
are the occupation of the system orbitals 1 and 2, respectively. They can have values or 1. The many-body level 
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scheme is sketched in Fig. 2. The full system density matrix has 16 components. In the reduced space, where FSC 
are eliminated, it has only six components. We order the vector given by the density matrix in the reduced space as 
P = (poo,oo , An, oi , Pio,io , Pii,n , Poi.io , Pio,oi)- Our QME, Eq. JTHJl, therefore reads 



with 



M = 



( -23?(/? n + foa) 
23t0 22 
29fy3n 


P21 + 

V 012 + 0*21 



23?a 22 
-2^{a 22 +0n) 


-a* 12 + 2 i 
-an + 02i 



p = Mp 



2Ka n 


-25R(an+/3 22 ) 
2$R/? 22 
-em + 0* 2 
-ali + /3i2 



(27) 





23?an 
23?a 22 
-25R(an + a 22 ) 

— Q!2l - a*2 
-a i2 - a 21 



«12 + a 21 
-a 21 + /3 i2 
-ai2 + /3|i 
-/3i2 - /?Ii 
-A" 




«21 - 
-«21 



'12 



l 12 



#2 

0*12 ~ 021 




(28) 



and X — + 0* x + a 22 + 22 s + i(ei — e 2 ). At steady state (p = 0) this equation can be transformed into a 4 x 4 
matrix equation for populations alone by including the effect of coherences into modified rates |36j . In the present 
work we do not consider the spin polarization of the electrons which has been used to study Pauli blockade |33j and 
magnetotransport |44| in QDs. Recently Gurvitz et al j3| have derived a QME to study the spin dependent coherence 
effects on electron transfer through a single QD. We notice that, in the high bias limit, our Eqs. 1)27(1 and ((28(1 are 
identical to their Eqs. (21) for the case of a single spin state (Sec. C in Ref. fSjp. 
The total charge of the system at time t is given by 



Q(t) = e <A|p(i)>= eTrNp(t) 
where N = c\c s is the number operator. The rate of change of the system charge is given by 

Q(t) = I L (t) + I R (t) , 
where II and Ir are the currents from the left and the right leads 

I x (t) = e <N\M X \p(t)>, X = L,R. 



(29) 



(30) 



(31) 



M.x is the contribution to the matrix M. from lead X so that M. — A4l + -Mr (Aix only contains terms in 128() 
corresponding to X lead). Since a(0) is related to the outfiux (influx) of electrons from (to) the system, we can 
further split the current as Ix = I™ + ^x" 



■IT*, where 



TOUt 

1 X 



™(t) = e^N\M x {0)W> 
(t) = e<^N\M x (a)\p{t)-> 



(32) 



Mx(a) [M x {0)] contains only those terms in (|2*S1) involving [0^ x ^]- At steady state (t — > oo), the currents 

from the left and from the right leads must be equal and opposite in sign, and the steady-state current is given by 

i s = i l = -i r . 

We solve Eq. 1(27(1 , by diagonalizing the matrix M. . 



|p»(t)-^C 1)e ^ t |r,», 



(33) 



where £ v (\r]^>) are the eigenvalues (eigenvectors) of M and C v =<C 7?|p(0) In the RWA [see Eq. (|2~T|) ] the 
off-diagonal terms a ss ' and SS ' are neglected and the population dynamics of Eq. 1(28(1 is then independent of the 
coherences and can be obtained analytically (see Appendix C). 

The steady state currents I m and I out are obtained from Eqs. ((32(1 and ((Cl(l 



o4(e s )a s 



= - 2e E 



ai{e s )0ss 



a^ s ( e s) + af s (e s ) 



h{e s ) 
0--h(e s )) , 



(34) 
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where a* s = nnx \T*(e s )\ 2 . Similar expressions can be obtained for the currents and l° R ut by interchanging L and 
R in Eq. ffiifr. Note that at steady state M L p = —M R p so that p = 0. Of course M L (a)p ^ —M L (f3)p. At steady 
state II = —Ir- We can therefore write for the steady state current I s = xlj, + — 1)Ir for arbitrary x. By choosing 
x = af s (e s ) I '(a,g S (e s ) + af s (e s )), I s can be written in a symmetric form [Tflll2lll3ll | 

'-*x £$g$ 3 v>M -/«(-»• <*» 

Since in the RWA the many-body density matrix is given by a product of single-orbital density matrices, the total 
steady state current is the sum of contributions from various orbitals. Note that Eci. (|35[) gives the steady state current 
within the RWA, which ignores the effects of coherences. In the model calculations presented in the next section, we 
shall discuss these results and the effects of coherences. 



V. NUMERICAL RESULTS 



We first solve Eq. (|27|l for the time-dependent density matrix in terms of the eigenvalues and eigenvectors of 
the matrix M. [see Eq. (|33H ]. We fix the system orbital energies, t\ — 5eV and 62 = 2eV and the temperature 
ksT = .2eV. In Fig. 3 we display the eigenvalue spectrum of M.. All eigenvalues have a real negative part 
representing an exponential decay. At long time, the system approaches the steady state corresponding to the zero 
eigenvalue. The two complex eigenvalues describe the two coherences. 

The time evolution of the populations [Eqs. [33] and coherences [Eqs. IC3| is shown in Figs. 4a and 4b, respectively. 
The two are decoupled in the RWA. Coherences show a damped oscillatory behavior and vanish at long times. 
The populations evolve exponentially and reach a steady state distribution described by the eigenvector with zero 
eigenvalue. 

For large bias (eV > £i) and identical left and right couplings (T S L = T s ), all many-body states have the same 
occupation. This is shown in Fig. 5. We assume that the chemical potential for the left lead increases with the bias 
while the right lead is held fixed at the Fermi energy po. When the bias is switched on, electrons start to move from 
the left to the right lead through the system. For po < and V = there arc no electrons in the system and the 
probability to find the system in the state |00) is one. This probability decreases as V is increased. Thus poo decreases 
with increasing bias. As the bias is scanned across higher energies, electrons start to fill the system. This gives rise 
to the step-wise change in the population in Fig. 5. As long as eV < ex, only states |00) and 1 01) are populated. For 
eV > ei, both system orbitals lie between po and po + eV, and all many-body states are equally populated. 

The steady state I/V characteristics of the system computed using Eqs. <|54l) and (|3*5|) are depicted in Fig. 6. The 
black solid curve shows the total current computed using Eq. (|35fl and dots (dash-dot) show the current 1™ (/£"*)• 
We note that /£"' is significant only at resonant energies where eV + po = e s . This can be explained as follows: when 
po + eV < e s , there are no electrons in the s'th orbital and hence = 0. For po + eV > e s , in order to move from 
s'th orbital to the lead, electrons need to work against the barrier generated by the chemical potential difference and 
hence the probability of back transfer is very small. At po + eV = e s , this barrier vanishes and electrons can move 
easily back to the left lead, giving rise to I out . 

We next compute the steady state charge on the molecule using Eq. I|C5|) . As the external bias is increased, 
different many-body states are populated and the charge on the system increases in steps, similar to the current. This 
is shown in Fig. 7 for different values of p . As the Fermi energy is increased the total system charge at steady state 
increases and the variation with the bias is decreased. Finally, when the Fermi energy is large enough so that all the 
many-body states are already populated at V = 0, the total charge (which is the maximum charge) on the molecule 
is 2e (both system orbitals are occupied) and is independent of the bias. 

We next study the effect of coherences by solving the QME (|27fl without invoking the RWA. In this case we need 
to diagonalize the full 6x6 matrix, A4 and we use Eq. (|31|l to compute currents numerically. In Figs. 8a and 
8b, we present the steady state currents with and without the coherences, respectively. The steady state coherences 
are shown in Fig. 8c. We note that due to coherences the backward current I out (dash-dot) does not vanish for 
eV + po ^ e s and is positive, although it is still maximum and negative at the resonances. This leads to the increase 
of the total current for smaller bias: coherences produce an effective potential that enhances the potential generated 
by the chemical potential difference between the leads. In Fig. 8d, we show the change in steady state currents due 
to coherences at different bias values. 
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VI. CONCLUSIONS 



The dynamics of a quantum system connected to two metal leads with different chemical potentials is calculated 
using projection operators which project the total many-body density matrix of the system into the system subspace 
corresponding to a fixed number of electrons n in Fock space. We derive a set of coupled dynamical equations for the 
n-dependent projected density matrix of the system. When summed over the different possible numbers of electrons 
n in the system we get a Redfield-like QME for the system many-body density matrix. Since we treat the leads as 
infinite electron reservoirs, coherences between the leads and the system are not possible. As a result, electron transfer 
between the leads and the system occurs in an incoherent way. This amounts to eliminating the coherences between 
system many-body states with different n (FSC) leading to a drastic reduction of the many-body system space. By 
studying the transient and steady state transport properties of a coupled QD system for arbitrary bias, we showed 
that coherences between the many body levels corresponding to a same n can affect the transport properties of the 
quantum system. In the limit of high bias our equation reduces to previously derived QMEs |28l l31| . By invoking the 
rotating wave approximation, we showed that the populations obey an independent birth and death master equation. 
In this limit, the QME can be solved analytically for an arbitrary number of orbitals since the system many-body 
density matrix is a direct product of the individual single-orbital density matrices. Both, in the high bias limit and in 
the RWA, our QME assumes the Lindblad form. We note that the full QME, Eq. i|10|). is not in the Lindblad form 
and may break positivity for strong lead-system couplings. Using the numerical solution of the QME for physically 
acceptable parameter range, we found that quantum coherences can modify the transport properties of the system. 
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APPENDIX A: DERIVATION OF THE QUANTUM MASTER EQUATION 

In order to compute the time dependence of p n (t) we start with the Liouville equation for the total density matrix, 
Pt- 

^ = -i[H T (t),p T (t)] (Al) 

where pr(t) represents the many-body density matrix and Ht is the coupling Hamiltonian, both in the interaction 
picture. 

p T {t) = e lHot p T (t)e- lHot (A2) 

where Hq = H s + Hl + Hr and pr(t) is in the Schrodinger picture evolving with the total Hamiltonian, H. The 
interaction picture operators are similarly defined by 

H T {t) = e lHot H T e- iHot . (A3) 

Substituting Ht from Eq. J5J in IjAlfl . multiplying by V n from both sides, taking a trace over the leads space and 
using Eq. we obtain the equation of motion for p n . 



dp n {t) 



= -iJ2 T ™ i A ™ (*) - B *v {t)]+h.c. ( A4) 



dt 
where 

A sv {t) = e ie '" tr Tj:i ead {c\c v Vn-ipT{t)V n } 

B su {t) = e^TTlead {VnPT(t)Vn+lctcu} , (A5) 

and e sl/ = e s — e u . In deriving Eq. ljA4fl . we have used the relations, 

V n c\c v = c\c v V n -i, V„clc s = 4c s "P n+1 . (A6) 
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Differentiating both sides of Eq. (|A5|) with respect to time and using Eq. O, we obtain 
dAsu{t) - ie^A^-ie^^e-^ 



dt 



{T S V [ctcs.Triead [c v cl,V n pr{t)Vn} - c\Tn ead \c v Vn-xpT{t)V n -xcl,} c, 



clTr /( 



(A7) 



dt 



It* , 



Cs'Tllead \ci<'P n+ lPT{t)'P n+ lC u } c\ - Tr; eQ(i {V n pT(t)Vncl,C^ C s >c\ 



(A8) 



This hierarchy involves successively higher coherences in Fock space. The first term in the r.h.s. of Eqs. I|A7(1 and 
(|A8I) represents the oscillatory time evolution due to the free molecule Hamiltonian. The other terms represent the 
evolution due to the coupling with the leads and involve populations and two-electron coherences in the molecule. 

We approximate each lead as a free electron gas described by the grand canonical density matrix psit) = PlPr, 
where p^ and pa are the density matrices for the left and the right leads with chemical potentials \il and /ir, 
respectively. We assume that the two leads have an infinite capacitance and are not affected by the weak coupling 
to the system. Both leads are therefore at equilibrium with their respective chemical potentials and the FSC in the 
leads vanish. This results in the loss of coherences between states with different number of electrons in the molecule 
and Eqs. (|A7f) and (|A8|) take the form 



dA sv (t) 
dt 

dB sv {t) 
dt 



jA sv (i) 



■~* e-^TX* [clc s ,p n (t)C^(t - t') - ctfr-WcV^it - t')] 



ie sl/ B sv (t) - ie 1 ^ e~ ie *'%*v> [c s > ~p n+1 {t)c\C v „, (t - t 1 ) - p n {t)c s> c\V^{t - t')] 



(A9) 



(A10) 



where C vu >(t — t') = Tr/ ead ^Cy{t)c\,(t')pB^ and V VV '{t — t') = Tr; ea( j jc£, (t')c v (t)pB\ are the correlation functions 

for the leads. 

The formal solution of Eqs. (|A9|) and (|A10(I is 



A av {t) 
B sv {t) 



-te 



■ f dt'yy^'T^, [cics.^it'^it-t^-cip^it'^v^it-t')] 

: / dty^e-^'T^, [c s ,p n+1 (t')c%„,(t-t') - p n (t')c s ,clV^(t-t')} 
Jo i i 



Since the leads are at equilibrium, their correlation functions are 

C vv ,(r) = <W(1- fv)z~ lt " T 
V vv ,{t) = 5 vv ,f„e~^ T 

where f v = [exp{/3(e„ — /!„)} + with \i v = /i£ or p^, v — I or r. 

Substituting Eqs. (|A11(1 and 1JA12I) in Eq. (IA4I) and making the change of variable, t — t' = r, we obtain, 

%T = ^ rdre^' i [a s y(r) Cs K-r)/'+ 1 (t-r)4-/3^(r)p"(t-T)^(-T)ct 

ss' ^ 

- a aa ,{T)clc s ,(-T)F(t-T)+0 as ,{T) C lr-\t-T)c s ,(-T)] + h.C. 

where e ss > = e s — e s > and we have used the notation 

<w(t) = ^r s ,T s %c^(r) = ^r s ,r;,(i-/, y ) e — ^ 

vv' v 



(All) 
(A12) 



(A13) 
(A14) 



(A15) 



(A16) 
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Transforming Eq. (|A1 5|l back to the Schrodinger picture, we get 

= -i[H 3 ,p n (t)]+Y / [ dT[a s AT)e- iHaT c 8 ,p n+1 {t-T)e iH ° T cl 

ss' 

- a s Ar)c\e-^ T c sl p n {t - r)e lH ° T ~ p ss , (r) e - iH ° T p n (t - r)c s ,e iH ^c\ 

+ 0„,(T)4e- iHoT P n - 1 (t-T)c a ,e iH ° T ]+h.c (A17) 

We next expand the density matrix p n (t — r) perturbatively in the coupling with the leads, 

p n (t-T) = e -iHa(t-r) p n e iH (t-r) +0 (T) 

= e iHoT p n (t)e- iHoT + 0(T) (A18) 

where 0(T) represents terms that depend on the leads-molecule coupling. Substituting Eq. (|A18J| in (|A17(I and 
keeping terms to second order in the coupling, we get, 

= -i[H s , p n (t)] + f dr [a ss , (r)<v (-r)p n+1 (t) c J - [5 SS , (r)p n (t)c s , {-r)c\ 

- a ss ,(T))4c sr (-T)p n (t)+P ss ,(T)cip n - 1 (t)c s 4-T)]+h.c. (A19) 

where c s (r) = e~ %e " T c s . Making the Markov approximation (assuming that the lead correlation time is short compared 
to the time evolution of p n ), the time integration in Eq. I|A19(1 can be extended to infinity and the equation becomes 
local in time. 

Substituting Eqs. I|A13|) in (|A19(1 and carrying out the time integration, we finally obtain Eq. I|10|) . Note that 
a similar derivation can be done in Liouville space [TH l45| by defining the Liouville space projection operator C n , 
C n pr = P n P B- 



APPENDIX B: QME IN THE LOCAL BASIS 



In this section, we recover Gurvitz's |2^| results starting from our QME I jlOjl for a QD system. Gurvitz considered 
a system of two QD connected in series between two leads. We denote the left and right QD by a and b respectively. 
We therefore need to transform the QME to the local basis representation. Let us define the unitary transformation 
matrix U which changes the system eigenbasis to local basis (denoted by indices where i — a,b and j — a, 6). We 

have / JUgjUja' = <W- The transformed Hamiltonian (JIJ in local basis then reads as 



with 



H = €i i c t c i + H £lC l Cl + £rC r Cr + \t iv c\c u + T* v clc 

ij I r iv 

i i ij 



Applying the unitary transformation, the QME can be transformed into the local basis set as 
dp n {t) 

ij ij 

- ocijc\c j p n {t) + I3 lj c\p n - 1 {t)c j + h.c. 



at 



where 



Pa = 



ss' 



(Bl) 



(B2) 



(B3) 



(B4) 
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Thus the QME structure remains the same. Note that even if we assume that the bath correlation function is diagonal 
in eigenstate, i.e. a ss i and (3 ss i are diagonal (which is equivalent to the rotating wave approximation, Sec. 1111(1 . so that 
the coherences become decoupled from the populations, in local basis however, since and /% are not diagonal [see 
Eq. (|B4() ]. the populations and coherences are still coupled. These coherences in the local basis, which are different 
from the coherences in eigenspace studied here, were analyzed by Gurvitz et al [28(. Our QME l|B3|) can be applied 
to Gurvitz's model of two QDs coupled in series, described by the Hamiltonian Q 



e ij c i c j 

^ eicjci , H R = t r c\c r 

I r 
^ ^r/(4cj + C v c\) 



(B5) 



where system Hamiltonian e aa — e a , e b b — e b and e a b = e ba = is the coupling between two dots. In the reduced 
Liouville space, where FSC are zero, we use the notation (n a nb\p\n' a Ji' b ) = p na n b ,n' a n' b m the local basis where n a (jib) is 
the occupation of QD a and b, respectively. The density matrix in the reduced Liouville space is given by the vector 
P= (/5oo,oo ! An, 01 ! Pio,io , Pii.n , Poi.io , /5io,oi)- Using Eq. (|B3|I . the time evolution of elements of the manv-bodv 
density matrix in the local basis is given by 



where 



M = 



( -23?(/3 aa + fi bb ) 
2U(3 bb 

23?/3aa 



Pba + P* ab 
\ Pab+PL 



2^a bb 

-2Vt{a bb + Paa) 



itt - a* ab + (3 ba 
-itt - a ab + (il a 



p = Mp, 



2Ma aa 


-23?(a aa + fab) 

-ifl - a ba + f3* ab 
i^o - a* ba + Pab 





23?a 

aa 

2Vta bb 
-25i(a Qa + a b b) 
-a ba - a* ab 
-dab - a* 



a a b + a ba 
iO - a* ba + (3 ab 
-iO - a ab + f3 ba 

~Pab ~ Pt a 

-X 




(B6) 



a ba + a* ab \ 

-ifl - a ba - fl* a b 
iQ, - a* ab + (3 ba 



—X* J 



(97) 



and X = a aa +f3* a +abb+f3bb + i{e a — £&)■ Note that Oiij — a^+af- and = P^+f3^. As done in Ref. |2S|. we assume 
a large external bias (/io + eV > e a , e b > po) so that the Fermi function for left lead is always 1 (occupied states) and 
that for the right lead is always zero (unoccupied states). In this limit electrons are always transferred from left to 
the right lead and the reverse transport vanishes, i.e. afj = f3fj = 0. We further assume that the bath correlation 



functions are real and that the lead's density of state is a constant (wide-band approximation). 
(1 A 16(1 in (|B4|I . it is easy to show that 



dij = 'l-ii n'l u I ' f, fa 
a a b= Pab =Pbb = 0, and T a = fL L , T b = fl R , we get 
a bb = 2irn R \n R \ 2 , f3 aa = 2Tm L \Q L \ 2 



2Trn L T z L T*f. 



Since = T b = 0, ct a 



Substituting Eqs. 



(B8) 



(B9) 



which is same as ^l(r) defined in Eq. (3.4) in Ref. [23]. The matrix jCi then simplifies to 



M = 



( 


-2K(3 aa 


25ia 6b 



















-2^(a bb +Paa) 










— iflg 




2^(3aa 








2®.a b b 


— iflg 









2^/3aa 





-2^a bb 








V 







— itto 





+«(e a - eb) - PL - a bb 








—i£lo 










-«(e a ^ eb) ~ Paa 



(BIO) 



l bb 
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We note that populations and coherences are coupled together only through the non-diagonal terms in the free 



(system) Hamiltonian £1 - The set of Eqs. (|B6[1 then can be written explicitly as 

^ = -2Kp aaP Ut) + ^a bhP ™+ l (t) (Bll) 

-gi = iQoipUt) - P21W) + 2m<aPoo 1 (t) + ma bb pl+ 1 {t) (B12) 

^ = iQ (p^(t)-pUt))-2^(a hb +/3 aa )p^(t) (B13) 

^ = 2^p aa p n 22 \t)-2^a bb p^{t) (B14) 

= ie ba pUt)+in (p n 11 (t)-pZ 2 (t))-2dt(a aa + f3 aa )pUt)- (B15) 



The only difference between Eas. ljBll() - (IB15|) and Eqs. (4.10a)-(4.10e) of Ref. is in the bookkeeping of electrons. 
In our case n is the number of electrons in the system whether in Ref. |28| n is the number of electron collected in 
the right lead. However, since we are in the large bias limit, the two are directly related and after summing over n, 
so that J2 n Pij = J2 n Pij 1 = Pv ' ^1 S ' ^11(1 - i|B15l) become identical to Gurvitz's equations. 



APPENDIX C: RWA SOLUTION FOR POPULATIONS AND COHERENCES 



In the RWA, populations and coherences evolve independently. As discussed in the main text, the population 
dynamics described by Eq. I|27[) obey a birth and death master equation. By diagonalizing the 4x4 generator of the 
population dynamics and using we get 

PWP11 P22 fill 

P11 Pn 

p WjW (t) = ^C 2 + — C 3 e- 2 ( Q11+/311 ^ - C 4 e- 2 ( Q22+& )* - C 1 e- 2 ( Qll+/3ll+Q22 + te ) t (CI) 
P22 P22 

where the coefficients C\ — C4 are related to the initial density matrix as follows 

011/322 



C 



2 



D 



C 3 = a " 22 (pn,n(0) + Pio,io(0)) - U D 22 (poo,oo(0) + Poi,oi(0)) 
C 4 = (poi,oi(0) + Pn,n(0)) - (Poo,oo(0) + Pio,io(0)) 

Ci = — {P11P22P00, oo(0) - an/?22Pio, io(0) - a 22 /3ii)poi, oi(0) + aiia 22 piiai(0)) , (C2) 

where p O o,oo(0) + Poi,oi(0) + Pio,io(0) + pii,n(0) = 1 and D = (a n + (3n)(a 22 + j3 22 ). 
The time-dependence of coherences is solely determined by the element X of matrix A4 

Poi.ioW = e— 12t e-( a "+^ +Q22 +' 322 )Voi,io(0). (C3) 

and pio.oi = Pqi w . When the bath correlation functions are real, coherences oscillate with Bohr frequency (£12) of 
the system. 

The steady state populations are given by 

Pll, 11 = J^PuP22, POO, 00 = -J^CillOi22 

1 1 1 1 

P01,01 = 7J«ll/?22, Pio.10 = -pAlQ!22- (C4) 
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Note that steady state coherences are zero. Using Eq. <|C4() in (|29|l it is then easy to show that the total steady state 
charge on the system is 
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Fig. 1: (a) Lead-system- lead configuration. a + (a~) and b + (b~) represent the charge transfer processes which 
change the number of electrons in the system due to interaction with the left (right) lead, (b) Energetics of the 
junction. //£ and hr are the chemical potentials of the left and right leads. E±, E 2 , £3 and E4 are the energies of 
the system many-body states. 

Fig. 2: The Four many-body states |ni,7i2) for a model system of two orbitals with occupations ni and 112 and 
energies ei and £2, respectively. N = Q, 1, 2 represents the total number of electrons in the system Fock space. Dashed 
and Solid lines denote the single-electron and many-electron states, respectively. E\,E%,E^ and E4 are the energies 
of the four many-body states. 

Fig. 3: The eigenvalue spectrum (in cV) of the matrix M for V = .1, Mo = 0, TL 1 = .01, TL 2 = .02, TR X = .03 
and TR2 = .04. All parameters are in units of eV. 

Fig. 4: (a) Time evolution of the populations (Eq. IC1JI for V = 2. Other parameters are same as in Fig. 2. Time 
is in units of h/eV. (b) Time evolution of the real (Repoi.io) and the imaginary parts (lmpoi,io) of coherences. 

Fig. 5: Steady state populations (Eq. IC4(1 for /io = 0. The left and right coupling are the same, TL\ = TR\ = 0.2 
and Tl2 — TR2 = 0.3. All parameters are in eV. 

Fig. 6: Current characteristics of the model system using Eqs. (|34f) and (|35|1 for parameter of Fig. 4. I m : dotted, 
jout. Jagged anc [ j s - solid curve. Current is in units of e 2 V/h. 

Fig. 7: Steady state charge (Q) on the system as a function of the applied bias (V) and the Fermi energy (^o)- 

Fig. 8: (a) Steady state currents obtained by solving Eqs. (|3lT) and 1(3*2^) and (b) without coherences, Eq. 
Comparing to (a), we note that in the absence of coherences, I° ut is significant only at the resonances and is always 
negative. (c)The steady state coherences in the system and (d) change in the steady state currents due to coherences. 
The couplings are TL X = .4, TL 2 = .7, TR X = .4, TL X = .2 and fj, = all in units of eV. 
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